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Abstract 

A  transf err  able,  polarizable  quantum  chemistry-based  force  field  has  been  developed  for 
hydrazinium  ,  monomethylhydrazinium,  and  dimethylhydrazinium  cations  in  combination  with 
the  nitrate,  azide,  dicyanamide,  and  5-azidotetrazolate  anions.  Quantum  chemistry  calculations 
were  performed  on  the  neutral  precursors,  ions  and  cation-anion  complexes  employing  aug-cc- 
pvDz  (cc-pvTz)  basis  functions  at  MP2  level  or  in  conjunction  with  M05-2X  density  functional. 
Inclusion  of  a  lone-pair  on  hydrazinium-based  cations  significantly  improved  ion  electrostatic 
description  and  prediction  of  the  crystal  structure  in  molecular  dynamics  (MD).  Seven  different 
ionic  systems  have  been  investigated:  [N2H5][N03],  [(CH3)N2H4] [N03],  [(CH3)2N2H3][N03], 
[N2H5][CN7],  [(CH3)N2H4][N3],  1(CH3)2N2H3][N3],  [N2H5][N(CN)2],  For  all  but  the 
[(CH3)2N2H3][N03]  and  [N2H5][N(CN)2],  QC  calculations  of  a  single,  gas  phase  ion  pair  predicts 
spontaneous  deprotonation  of  the  cation. 

Crystal  lattice  parameters  obtained  from  MD  simulations  were  compared  with  experiments  for 
ionic  crystals  of  these  seven  systems,  with  the  experimentally  determined  crystal  structure  of 
[N2H5][N(CN)2]  also  presented  here,  enabling  comparison  of  simulation  and  experiment  for  that 
compound.  In  general,  MD  simulations  predicted  crystal  lattice  vectors/angles  (volumes)  within 
a  5%  (3%)  absolute  margin  of  error  from  experiments,  with  outlying  values  of  5-6.6  %  for  three 
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crystals  [(CH3)N2H4][N3],  [N2H5][N03],  [(CH3)N2H4][N03]  with  combinations  of  particularly 
small  anions  and/or  cations.  Structural  comparisons  between  ionic  materials  in  the  liquid  and 
crystalline  state  are  made,  including  the  observation  of  two  crystalline  systems  where  the 
crystalline  state  induces  conformational  changes  in  the  methylated  hydrazinium  cations  compared 
to  the  gas  phase  and  liquid  states.  The  matrices  of  elastic  constants  were  extracted  from  MD 
simulations  for  all  seven  ionic  crystals  and  correlated  with  the  structural  motifs  of  ion  interactions 
in  the  crystals. 

Introduction 

Hydrazine,  mixtures  of  hydrazine,  or  derivatives  of  hydrazine  have  been  utilized  for  wide- 
ranging  applications  including  polymer  synthesis  and  processing,  metal  recovery  and  removal, 
pharmaceutical  and  agricultural  applications,  and  bio-active  applications  such  as  herbicides  and 
pesticides.1  Perhaps  the  highest  profile  use  for  the  hydrazine-based  compounds,  however,  is 
provided  by  its  longstanding  role  as  a  primary  constituent  of  fuel  systems  for  aerospace 
applications.  Much  research  has  been  performed  to  elucidate  the  behavior  of  hydrazine  and  its 
methylated  derivatives,  monomethylhydrazine  (MMH)  and  unsymmetric  dimethylhydrazine 
(UDMH)  by  NASA2'3  and  others  in  order  to  understand  and  control  the  conversion  of  these  high 
energy  chemicals  into  utilizable  energy  sources  for  aerospace  applications.  Unfortunately,  the 
properties  that  make  hydrazine  and  derivatives  desirable  compounds  for  these  applications 
(namely  high  energy  density,  low  reaction  threshold,  and  high  volatility  associated  with  the  liquid 
state)  also  contribute  to  its  potential  for  harm  during  uncontrolled  release.  Hydrazine,  MMH,  and 
UDMH  have  been  shown  to  be  carcinogenic4'5  in  laboratory  tests  involving  rodents,  and  are 
suspected  to  have  similar  effects  on  humans  via  oral  or  vapor  exposure. 

Recently,  ionic  liquids  (ILs)  have  attracted  significant  attention  as  solvent  replacements6'7  and 
alternatives  to  the  hydrazine-based  hypergolic  propellant.8 11  Due  to  their  ionic  nature,  IL  vapor 


Distribution  A:  Approved  for  public  release;  distribution  unlimited 


pressure  is  low,  resulting  in  a  reduced  environmental  impact.  Moreover,  existence  of  numerous 
cation-anion  combinations  provides  an  avenue  for  optimization  of  the  thermodynamic  properties 
and  reaction  pathways.  Experimental  synthesis  and  characterization  of  novel  high  nitrogen 
containing  ionic  liquids  and  crystals  have  been  hindered  by  reactivity,  adding  to  the  challenges 
and  safety  risk  for  experimental  investigation.  In  cases  such  as  these,  the  ability  to  accurately 
predict/describe  thermodynamic  response  of  these  compounds  without  the  associated  risks  of 
experimental  measurements  is  of  obvious  utility.  While  there  are  many  important  properties  that 
simulations  have  difficulty  addressing,  thermodynamic  behavior  of  the  material  can  be 
quantitatively  described  given  an  accurate  force  field  description  and  proper  simulation 
methodology.  We  have  recently  demonstrated  the  ability  of  our  force  field  development  and 
simulation  process  to  provide  simulation  models  that  can  accurately  reproduce  and  predict 
density,  heat  of  vaporization,  ion  diffusion  coefficient,  conductivity  and  viscosity  of  a  diverse 
class  of  ionic  liquids.12  The  same  force  field  accurately  predicted  ionic  crystal  lattice  parameters 
for  four  ionic  crystals12  as  well  as  prediction  of  lattice  parameters,  conductivity,  and  melting  point 
of  a  plastic  crystalline  system.13  Recently,  however,  Maginn’s  group  has  shown14  that  it  is  quite 
difficult  to  achieve  an  accurate  and  consistent  description  of  the  modified  hydrazinium  nitrate 
ionic  crystals.  In  this  work,  we  undertake  the  challenge  of  developing  a  systematic  force  field  that 
will  describe  crystal  structures  for  these  modified  hydrazinium  nitrate  ionic  materials. 
Furthermore,  the  developed  polarizable  force  field  will  be  extended  to  other  anions  such  as  azide 
(N3"),  dicyanamide  (N(CN)2‘ )  and  5-Azidotetrazolate  (CN7')  in  addition  to  NO3  and  will  be 
validated  using  seven  high  nitrogen  content  ionic  crystals.  After  force  field  validation,  we  use 
MD  simulations  employing  the  developed  force  field  to  perform  initial  investigations  of  the 
mechanical  properties  in  the  crystalline  system  and  investigate  the  structural  interactions  of 
cations  and  anions  in  the  crystalline  (and  in  some  cases  liquid)  state,  focusing  on  the  differences 
of  ion  packing  in  the  crystalline  and  liquid  states. 
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The  primary  difficulty  of  developing  a  force  field  for  the  high  nitrogen  content  ionic  liquids 
and  crystals  is  a  lack  of  thermodynamic  information  in  the  liquid  phase,  which  is  often  unstable 
for  hydrazinium-based  ionic  liquids.  Therefore,  we  rely  heavily  on  the  crystal  structures  of  the 
hydrazinium-based  ionic  crystals  for  force  field  validation.  Starting  from  a  description  of  the 
neutral  hydrazines  as  a  baseline,  we  refine  our  force  field  models  to  produce  the  best  overall 
description  over  the  entirety  of  the  crystal  structures  utilizing  a  single,  consistent  force  field,  with 
the  idea  that  utilizing  a  consistent  set  of  potential  interactions  for  a  number  of  ionic  crystals  with 
different  crystal  symmetries,  melting  temperatures,  and  anions  will  provide  sufficient  statistical 
variance  over  the  specific  sampled  interactions  of  the  hydrazinium-based  ions  to  allow  a  good 
parameterization  of  a  hydrazinium  force  field  for  use  in  the  generic  case. 

To  this  end,  we  have  investigated  three  substituted  forms  of  hydrazine,  in  both  their  neutral 
and  singly  charged  cationic  states:  Hydrazine  and  its  cation  (hydrazinium,  N2H5+), 
monomethylhydrazine  (MMH)  and  its  cation  (methyl  hydrazinium,  (CH3)N2H4+),  and  the 
unsynmietric  dimethylhydrazine  (UDMH)  and  it’s  cation  (1,1-N-Dimethylhydrazinium, 
(CH3)2N2H3+).  Representations  of  the  individual  neutral  molecules  are  presented  as  Figure  1A, 
while  the  cationic  forms  are  presented  as  Figure  IB.  In  the  case  of  the  cationic  forms,  we  have 
investigated  a  number  of  ionic  crystals  which  utilize  four  common,  singly  charged  nitrogen 
containing  anions:  Nitrate  (N03‘),  azide  (N3  ),  dicyanamide  (N(CN)2‘)  and  5-Azidotetrazolate 
(CN7'),  shown  in  Figure  1C.  Note,  that  the  CN7  anion  has  one  of  the  highest  known  nitrogen 
content  and  was  only  recently  characterized  in  detail.13  For  the  NO?  anion,  crystals  with  all  three 
cationic  hydrazinium  forms  were  available  and  will  be  utilized  for  the  force  field  validation, 
while  both  the  N(CN)2‘  and  CN7  anion  were  investigated  with  only  the  N2H5+  cation.  For  NY, 
only  crystal  structures  involving  the  methylated  cations  ((CH3)N2H4+  and  (CH3)2N2H3+)  were 
available.16 17  For  the  [N2H5][N(CN)2]  system,  X-Ray  diffraction  analysis  has  been  performed 
determine  the  crystal  structure  and  allow  for  simulation  comparison.  The  pertinent  experimental 
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details  for  the  determination  of  the  crystal  structure  are  presented  in  the  supplemental  information 
for  this  paper,  while  the  crystal  structure  itself  is  available  through  the  Cambridge 
Crystallographic  Data  Center  (CCDC)  as  structure  XXXXXX. 

The  remainder  of  the  paper  is  organized  as  follows:  In  Section  II  the  force  field  development 
is  outlined,  while  in  Section  III  we  discuss  simulation  methodology.  Section  IV  presents  the 
results,  while  interpretation  and  directions  for  future  work  are  provided  in  section  V. 

II.  Force  Field  Development 

Force  Field  Functional  Form.  The  following  form  of  the  force  field  relating  the  potential 
energy  U™  (r)  to  atomic  coordinates  r  for  the  ensemble  of  atoms  has  been  chosen.  It  is  split 
into  nonbonded  Unb{t)  and  bonded  contributions  as  follows: 

UT“(r)  =  U"(r)+'2uM(eill)+  J  <» 

Bends  Dihedrals  Improper 

Dihedrals 


where  the  sums  are  over  all  bends,  dihedrals,  and  improper  dihedrals  in  the  system.  The 
contributions  of  each  of  these  terms,  respectively,  are  calculated  as: 


T  jBend  ( t  \  \  ^  /  Bend  t  r\ 

U  (eijk)  =  -kaflY  (eijk 


(2) 


U 


Dihedral 


(3) 

(4) 


where  0i/k  and  0f/k  are  the  instantaneous  and  natural  bending  angles  for  atoms  i,j,  and  k\  (f)ljkl 
is  the  dihedral  angle  for  atoms  i,j,  k,  and  /;  and  r//"^  is  the  out-of-plane  bending  angle  for  an  sp2 
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center  at  atom  j.  The  strength  of  these  interactions  is  characterized  by  the  corresponding  force 
constants:  kB^d  ,  k^e6d'al , and  k^d  respectively,  where  the  subscripts  a,  />,  y,  <5  represent  the 

atom  type  for  atoms  i,j,  k,  and  1  respectively.  There  are  no  bond  energy  terms  due  to  the  fact  that 
bonds  are  always  constrained  at  a  bond  specific  fixed  distance. 

The  nonbonded  energy  [{/^(r)]  is  generated  by  summing  the  two-body  repulsion  and 
dispersion  energy  terms  [  £/SD(l")],  the  energy  due  to  interactions  of  fixed  charges  [  Utoul (t)], 
and  the  polarization  energy  [  UIoI{t)],  arising  from  the  interaction  between  induced  dipoles  with 
both  fixed  charges  and  other  induced  dipoles: 


UNB  (r)  =  URD  (r)  +  UCoul  (r)  +  UPo'(r) 


where  /i(  =  is  an  induced  dipole  at  force  center  i,  a.  is  the  isotropic  atomic 

polarizability,  Ej°'  is  the  total  electrostatic  field  at  the  atomic  site  i  due  to  permanent  charges  qj 
and  induced  dipoles  Jlj,  f0  is  the  dielectric  permittivity  of  vacuum,  £(°  is  the  electric  field  due  to 
fixed  charges  only,  Aap  and  Ba/i  are  the  repulsion  parameters  and  Ca/}  the  dispersion  parameter 

for  interaction  between  atoms  i  and  j  with  atom  types  a  and  [>.  The  term  D^l  2/ r.  j  ,  with 

D  —  5  x  10“5  kcal/mol  for  all  pair  interactions  not  involving  a  lone  pair  force  center,  is 
essentially  zero  at  typical  nonbonded  atomic  separations,  but  becomes  the  dominant  term  at  ry  <  1 
A,  ensuring  that  URD( I")  is  repulsive  at  distances  much  smaller  than  the  size  of  an  atom. 
Intramolecular  nonbonded  interactions  are  included  for  atoms  separated  by  three  or  more  covalent 
bonds.  We  have  utilized  Thole  screening18  ( aT  =  0.2),  which  locally  “smears”  induced  dipoles 
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over  a  short  spatial  range,  in  order  to  prevent  the  so-called  “polarization  catastrophe”  from 
occurring.  Additionally,  a  reduction  by  0.8  is  applied  to  intramolecular  interactions  between  an 
induced  dipole  and  a  partial  charge  separated  by  three  bonds.  Finally,  for  heteroatom 
interactions,  the  modified  Waldman-Hagler  combining  rules  were  used: 


, -  B6 

Aij  =  V AaAii  />/>< 

a  jj 

(  2  r 

Bij  =  n-6  ,  n-6  ^ 

\Bii  +BJj) 

c ..  =  Jc  c 

v  V  "  jj 

These  combining  rules  have  been  successfully  used  by  us  for  simulations  of  liquids,  polymers, 
electrolytes,  and  ionic  liquids12.  A,  B,  and  C  parameters  can  be  expressed  in  terms  of  potential 
well  depth  e,  the  interatomic  separation  at  the  minimum  R  ,  and  the  steepness  parameter  A  as 
given  by  eq.  7-9: 


A  =6£(expA)/(A-6) 

(V) 

B  =  A /R* 

(8) 

C  =  eX(r*)6 /{A -6) 

(9) 

Partial  Atomic  Charges.  Partial  atomic  charges  are  based  on  the  electrostatic  potentials  of  the 
gas  phase  molecules  calculated  on  a  grid  of  evenly  spaced  points  (~105  points  per  molecule) 
around  the  lowest  energy  gas-phase  conformations,  with  a  small  weighting  also  given  to  the 
dipole  moment  ( /7;)  of  the  molecule.  Cationic  and  neutral  molecular  species  were  based  on 
quantum  results  calculated  at  the  MP2/aug-cc-pvDz  level,  while  anionic  species  were  calculated 
at  the  MP2/cc-pvTz  level,  with  charges  assigned  in  the  QC  results  via  the  CHelpG  charge 
distribution  methodology.19  Charge-bond  increments  are  used  to  calculate  partial  atomic  charges, 
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where  the  value  ( qt)  of  the  partial  charge  positioned  on  atom  i  is  calculated  as  the  sum  of  all 
charge-bond  increments  that  involve  atom  i: 


N Bond 

V!=28'J  d°) 

7=1 

with  dy  the  charge-bond  increment  connecting  an  atom  of  type  i  to  an  atom  of  type  j.  The  set 

of  charge-bond  increments  for  the  overall  molecule  is  determined  by  minimizing  the  objective 
function: 


M 

1=1 

where  (pff  and  (j)1,  '  are  the  electrostatic  potential  for  the  /th  molecule  (or  complex)  at  the  jth 

grid  point  from  QC  and  the  developed  force  field  (FF)  respectively,  Jlf1  and  J1FF  arc  the  dipole 

moments,  and  of  and  of  are  the  relative  weights  for  fitting  the  electrostatic  potential  and  the 
dipole  moments  respectively.  The  electrostatic  potential  weighting  was  set  to  1.0,  with  the  dipole 
weighting  set  to  0.1  for  all  systems.  All  points  within  1.8  A  of  an  oxygen  atom,  1.5  A  of  a 
hydrogen  atom,  2.5  A  of  a  carbon  atom,  or  2  A  of  a  nitrogen  atom  or  further  from  any  atom  than 
4  A  were  excluded  from  the  fitting.  For  the  dicyanamide,  nitrate,  and  azide  anions,  this  process  is 
sufficient  to  yield  a  good  quality  approximation  of  the  electrostatic  field  due  to  the  molecules. 
However,  for  all  the  hydrazinium  cations  and  the  5-Azidotetrazolate,  fitting  charge  increments 
only  to  the  explicit  atoms  in  the  molecule  yielded  unsatisfactory  results  in  terms  of  accurate 
reproduction  of  the  molecular  electrostatic  field.  In  both  cases,  the  solution  was  to  add  explicit 
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off-atom  force  centers  to  act  as  additional  degrees  of  freedom  for  the  charge  fitting  routine  in 
order  to  improve  the  electrostatic  potential  fit.  We  therefore  add  dummy  force  centers  to 
represent  “lone  pairs”,  with  the  general  approach  utilized  explained  in  detail  elsewhere.20,21  The 
extended  charge  center  is  non-polarizable  and  has  zero  mass,  contributing  only  to  electrostatic 
non-bonded  interactions  via  acting  as  a  static  partial  charge.  The  electrostatic  force  induced  on 
the  dummy  force  center  is  distributed  to  its  attached  tri-atomic  basis  using  the  chain  rule.22 

The  importance  of  including  the  lone  pairs  for  accurate  reproduction  of  the  local  electrostatic 
field  of  the  hydrazinium  cations  is  demonstrated  in  Figure  2.  When  the  explicit  lone-pairs  are 
foregone,  the  potential  near  the  NFF  group  opposite  the  positively  charged  nitrogen  is  smooth,  as 
shown  in  Fig.  2a.  The  quantum  chemistry,  on  the  other  hand,  shows  a  marked  depression  in  local 
field  (which  shows  up  as  a  “well”  in  a  constant  potential  isosurface  as  plotted  in  figure  2)  in  a 
position  which  is  commensurate  with  the  location  one  would  expect  from  a  lone  pair  in  a 
tetrahedral  configuration  attached  to  the  nitrogen  in  the  NFF  group.  In  Fig.  2b,  we  present  a 
comparison  of  the  electrostatic  potential  surfaces  for  our  hydrazinium  cation  with  the  addition  of 
two  dummy  force  centers.  With  the  addition  of  the  lone  pair  force  centers,  the  simulation 
electrostatic  field  becomes  far  better  at  approximating  the  quantum  chemistry  potential, 
approximating  the  potential  “well”  at  a  constant  electrostatic  potential  with  the  same  value  as 
depicted  in  Fig.  2a  quite  well.  It  should  be  noted  that  the  protruding  “finger”  of  quantum 
potential  isosurface  that  extends  further  toward  the  nitrogen  in  Fig.  2b  is  due  to  the  fact  that  the 
range  of  the  QC  surface  actually  impinges  on  the  excluded  radius  (2  A)  of  the  nitrogen  atom  in 
the  partial  charge  calculation  (placing  it  well  into  the  repulsive  core  of  the  nitrogen  atom).  In  the 
force  field,  the  grid  points  closer  than  2.0  A  to  the  nitrogen  were  excluded  from  the  fit  and  are 
therefore  not  shown  in  Figure  2b. 

Introduction  of  the  lone  pairs  to  the  NH2  group  has  a  substantial  effect  on  the  overall  quality 
of  the  electrostatic  potential  fit  in  quantitative  terms,  with  the  best  fit  obtained  without  lone  pairs 
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yielding  approximately  1.49,  1.23,  and  0.98  kcal/mol  mean  square  deviation  observed  for  N2H5+, 
(CH3)N2H4+,  and  (CH3)2N2H3+.  With  the  addition  of  both  lone  pair  force  centers,  those  values 
reduce  to  0.77,  0.77,  and  0.78  kcal/mol,  respectively.  Interestingly,  if  we  simply  add  the  single 
(outward)  lone  pair  in  the  direction  of  the  observed  QC  “well”  and  calculate  the  charge  fit,  the 
results  are  only  marginally  better  than  for  no  lone  pairs.  This  resulting  necessity  of  an  internally 
charged  “phantom”  lone  pair  is  similar  to  the  significantly  improved  charge  fit  which  occurs  with 
a  non-intuitively  placed  extended  force  center  in  the  TIP4P  water  model.23 

In  the  5-Azidotetrazolate  anion  as  well,  lone  pairs  significantly  improve  the  reproduction  of 
the  electrostatic  field  near  the  molecule.  Here,  however,  the  difference  is  subtler  in  nature.  Initial 
predictions  of  the  electrostatic  field  around  the  tetrazolate  ring  without  the  lone  pair  force  centers 
yields  a  field  which  is  too  round  in  character,  as  demonstrated  by  the  red  surface  in  Figure  3.  Due 
to  the  aromatic  nature  of  the  nitrogens  in  the  ring  (as  well  as  the  resonant  nature  of  the  azide 
nitrogen  connected  to  the  ring's  single  carbon  atom)  we  have  added  planar  lone  pairs  to  all  of  the 
ring  nitrogens  as  well  as  the  azide  nitrogen  connected  to  the  ring  carbon.  The  resultant 
electrostatic  fit  to  the  quantum  chemistry  is  shown  by  the  blue  surface  in  Fig.  3,  and  can  be  seen 
to  be  generally  excellent,  with  only  a  small  bit  of  the  “cusp”  near  the  azide-carbon  bond  not  well 
accounted  for.  For  this  system,  the  initial  mean  square  deviation  between  force  field  and  QC 
calculations  is  3.14  kcal/mol,  while  the  inclusion  of  the  explicit  lone  pairs  drops  this  value  to  1.00 
kcal/mol. 

We  did  not  fit  charges  for  the  cation-anion  pairs  because  proton  transfer  was  found  via  QC 
calculation  at  the  MP2/aug-cc-pvDz  level  of  theory  for  all  ion  pairs  except  the 
[(CH3)2N2H3]  [N03]  and  [N2H5][N(CN)2]  systems.  In  the  former  exception  the  nitrate  anion  and 
dimethylhydrazinium  cation  jointly  share  the  proton  in  the  minimum  energy  configuration,  while 
only  in  the  hydrazinium/dicyanamide  system  do  the  cation  and  anion  retain  their  ionized  nature  in 


Distribution  A:  Approved  for  public  release;  distribution  unlimited 


the  gas  phase.  For  the  same  reason,  we  did  not  validate  ability  of  the  force  field  to  reproduce  the 
cation-anion  binding  energy  obtained  from  quantum  chemistry  calculations. 

For  the  [N2H5][N03],  we  have  run  further  calculations  using  the  polarizable  continuum  model, 
calculating  [N2H5][N03]  pair  interactions  within  a  cavity  which  is  in  turn  embedded  in  a 
continuum  model  background  with  the  dielectric  properties  of  acetone.  In  the  cavity  within 
acetone,  a  hydrazine/nitric  acid  pair  at  the  interaction  distances  and  geometries  generated  from 
the  isolated  gas  phase  QC  calculation  undergoes  ionization,  regenerating  the  expected 
[N2H5][N03]  ion  pairing. 

Atomic  Polarizabilities.  Isotropic  atomic  polarizabilities  are  calculated  from  both  the  gas 
phase  molecular  dipole  moment  and  explicit  binding  energy  QC  calculations  between  the 
molecule  of  interest  and  a  test  charge  performed  at  the  MP2/aug-cc-pvDZ  (MP2/cc-pvTz)  level  of 
theory  for  neutral  and  cationic  (anionic)  molecules.  Due  to  the  transferable  nature  of  the  force 
field,  it  is  generally  only  necessary  to  calculate  the  polarizability  of  a  specific  atom  type  once  for 
each  general  chemical  environment  it  occupies,  after  which  this  polarizability  is  used  transferably 
for  other  atoms  in  similar  chemical  environments.  This  approach  has  yielded  good  results  for  a 
number  of  previous  investigations.12'24  25  Polarization  interactions  between  atoms  separated  by 
one  or  two  bonds  are  assumed  to  be  subsumed  within  the  bond  and  bend  respectively,  and  are 
therefore  excluded  from  our  calculations. 

In  this  work  the  practical  upshot  of  this  approach  is  that  we  have  calculated  the  polarizability 
of  the  nitrogen  in  the  neutral  hydrazine,  and  transferred  the  atomic  polarizability  to  the  cationic 
hydrazinium  nitrogen  atoms.  In  combination  with  the  charge  fitting  outlined  previously,  this 
yielded  predictions  for  molecular  polarizability  of  the  cations  which  were  within  5%  of  the  QC 
values,  with  N2H5+  having  a  value  of  2.48  A3  (QC  =  2.53  A3),  (CH3)N2H4+  at  4.13  A3  (QC  =  4.35 
A3),  and  (CH3)2N2H3+  showing  5.80  A3  from  force  field  calculations  vs.  a  QC  value  of  6.00. 
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For  the  anionic  species,  we  fit  only  the  azide  and  CN7  anions,  since  the  N(CN)2"  and  NO3' 
anions  were  previously  fit.12  We  first  fit  the  polarizability  of  the  azide  anion  using  a  teardrop 
shaped  charge  path  both  along  and  at  a  45  degree  angle  to  the  principle  axis  of  the  molecule, 
yielding  a  calculated  molecular  polarizability  of  3.37  A3  for  the  force  field  vs.  3.39  A3  for  the  QC. 
We  then  utilized  these  polarizabilities  directly  for  the  azide  tail  of  the  CN7  molecule,  as  well  as 
transferring  the  terminal  nitrogen  polarizability  from  the  azide  to  the  CN7  ring  nitrogens  since 
both  molecules  are  comprised  of  resonant  anionic  nitrogen  atoms.  The  resultant  molecular 
polarizabilities  for  the  CN7  ion  predicted  by  the  force  field  was  8.65  A3  compared  to  the  QC 
value  of  11.35  A3.  Attempts  to  increase  the  molecular  polarizability  to  more  accurately  represent 
the  quantum  values  provided  a  dilemma:  We  could  easily  increase  the  value  of  the  polarizability, 
but  only  at  the  expense  of  the  quality  of  fit  of  the  overall  molecular  dipole.  The  as-transferred 
parameters  were  found  to  yield  the  best  fit  to  the  dipole  of  the  overall  molecule,  yielding  under¬ 
prediction  of  the  molecular  polarizability  by  ~25%.  This  level  of  reduction  in  the  overall 
molecular  polarizability  has  been  found  to  be  necessary  in  previous  simulations  that  were  able  to 
accurately  reproduce  the  crystal  structure  of  the  TATB  molecular  crystal.26  It  is  noteworthy  that 
TATB  also  represents  a  molecular  system  comprised  of  resonant  ring  structures  with  the  presence 
of  strong  hydrogen  bonding. 

Valence  and  Dihedral  Parameters.  Bond  lengths  were  constrained  for  the  MD  simulations, 
with  values  fixed  to  reproduce  those  in  single-molecule  QC  calculations  at  the  MP2/aug-cc-pvDz 
(neutral  and  cations)  and  MP2/cc-pvTz  (anions)  theory  level.  The  QC  calculations  using  triple-^ 
basis  sets  for  anions  were  chosen  over  the  aug-cc-pvDz  basis  set  due  to  their  ability  to 
consistently  better  represent  the  bond  lengths  observed  via  X-ray  measurement  of  the  ionic 
crystals.  The  natural  bending  angles  ( d-jk )  were  fit  to  the  optimized  QC  geometries,  with  angles 

in  common  chemical  environments  constrained  to  have  the  same  force  field  parameters  across 

instances  where  the  angle’s  central  atom  is  bonded  to  the  same  elements,  in  the  same  nominal 
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charge  state,  in  the  same  or  symmetrically  equivalent  orders  (i.e.  one  set  of  common  angular 
parameters  for  the  H-N-H  angle  of  the  NH2  group  bonded  to  the  “neutral”  nitrogen  in  the 
hydrazines  regardless  of  methylation  state,  another  set  for  the  cationic  hydraziniums  where  the 
NH2  group  is  bonded  to  the  formally  +1  charged  nitrogen,  with  parameters  again  equal  across  all 
three  methylated  derivatives). 

The  torsional  (dihedral)  parameters  are  fit  to  energies  from  MP2/aug-cc-pvDz//MP2/cc-pvTz 
level  scans  of  the  dihedral  angle(s)  fixed  in  ten  degree  increments  up  to  the  maximum  symmetry- 
reduced  angle,  with  the  remainder  of  the  molecular  degrees  of  freedom  allowed  to  relax.  In 
general,  since  the  torsional  description  of  the  chain  is  intimately  linked  to  the  description  of  the 
individual  bends,  we  pursued  a  self-consistent  state  for  the  parameterization  of  both  the  bending 
and  dihedral  angles,  with  fitting  alternating  between  first  one  and  then  the  other  sets  of 
parameters  until  the  parameters  became  stable  with  respect  to  fitting  the  complementary  set  of 
parameters.  We  found  that  achieving  this  level  of  stability  in  the  description  of  the  bending  and 
torsional  characterization  led  not  only  to  geometrically  consistent  descriptions  of  the  isolated 
molecule,  but  also  to  a  valence  description  in  which  the  out-of-plane  (improper)  torsional 
contribution  was  minimal,  with  the  related  energetic  constant  essentially  remaining  zero.  This 
latter  is  especially  important  since  the  NH2  group  is  capable  of  undergoing  inversion  at  high 
temperature,  which  circumvents  the  full  torsional  energy  in  certain  systems,  and  it  is  therefore 
important  to  accurately  capture  the  inversion  profile  of  the  improper  torsion  to  ensure  proper 
energetics  of  the  hydrazines  and  hydraziniums. 

Repulsion-Dispersion  Parameters.  Fitting  repulsion-dispersion  (R/D)  interaction  parameters 
(A,B,C  in  eq.  6  or  e,  X,  R  in  eqs.  7-9)  is  the  most  challenging  and  controversial  part  of  force 
field  development  since  the  parameters  can  be  obtained  by  fitting  to  a  number  of  different  data 
sets  yielding  somewhat  different  potentials.  The  most  often  used  data  for  fitting  the  repulsion 
dispersion  parameters  are  crystal  structures  together  with  sublimation  energies,  liquid  densities, 
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and  heats  of  vaporization  (Hvap),  vapor-liquid  equilibrium  data,  and  gas-phase  dimer  binding 
energies  obtained  from  QC  calculations.  In  previous  work  we  have  fit  the  R/D  parameters  to  the 
density  and  Hvap  of  nonionic  liquids,  while  the  density  and  dynamic  properties  were  used  for  IL 
specific  parameters.12  While  this  approach  has  proven  to  work  very  well,  it  is  intractable  for  our 
current  investigation  because  of  the  dearth  of  information  concerning  dynamic  properties  of  most 
of  the  ILs  investigated.  In  general,  due  to  the  energetic  nature  of  these  compounds,  little 
experimental  data  is  known  regarding  dynamic  properties  for  the  ionic  forms  of  the  compounds  of 
interest  in  the  liquid  state.  Thus,  we  have  fit  R/D  potentials  of  the  ionic  forms  via  transferring 
them  from  the  non-ionic  precursors,  which  have  in  turn  been  fit  primarily  based  on  reproducing 
the  liquid  state  densities  and  Hvap  for  clear  precursors  where  available  (e.g.  the  substituted 
hydrazines)  or  for  an  unambiguously  functionalized  derivative  (e.g.  1-Azidobutane  for  the  azide). 
The  sole  exception  to  this  was  the  fitting  of  the  (heavily  screened)  N+  homoatomic  potential, 
which  has  been  generated  by  rescaling  the  neutral  nitrogen  R/D  potential  as  discussed 
previously.12  For  DCA  and  NO3 ,  full  parameterizations  of  good  accuracy  were  available  from 
previous  efforts  with  differing  cationic  species.12'13 

TIT.  Simulation  methodology 

Simulations  in  both  the  crystalline  and  liquid  phases  have  been  carried  out  utilizing  the 

Lucretius27  molecular  simulation  package.  Periodic  boundary  conditions  were  used  in  all 

simulations,  with  covalent  bond  lengths  constrained  using  the  velocity- Verlet  form  of  the 

SHAKE  algorithm.28  The  Ewald  summation  method  was  used  for  the  treatment  of  long  range 

electrostatic  forces  between  partial  charges,  as  well  as  between  partial  charges  and  induced 

atomic  isotropic  dipoles,  while  direct  induced  dipole-induced  dipole  interactions  were  calculated 

only  via  an  iterative  realspace  summation.  The  number  of  reciprocal  Ewald  vectors  was  chosen 

to  satisfy  the  condition^  >  Ljljl  with  k/  the  number  of  Ewald  vectors  and  L.  the  total  length 

of  the  simulation  cell  along  lattice  vector  i,  with  a=0.23  used  for  the  real-space  decay  factor.  The 
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number  of  ion  pairs  varied  from  160  to  384,  yielding  periodic  cell  lengths  of  30-40  A,  with  the 
particular  value  between  these  limits  dictated  by  the  size  of  the  crystalline  unit  cell  in  most  cases. 
The  MD  simulations  employed  a  multi-timestep  reversible  reference  system  propagator 
algorithm29,  with  a  time  step  of  0.5fs  for  bonds  and  bend  vibrations  and  a  time  step  of  1.0  fs  for 
dihedral,  out-of-plane  deformation,  repulsion-dispersion  (van  der  Waals)  interactions  and  the  real 
part  of  Ewald  within  7  A  cutoff.  A  timestep  of  2.0  fs  was  used  for  all  interactions  within  a  1 1  A 
cutoff  including  the  reciprocal  part  of  Ewald  and  many-body  polarization.  A  fourth  order 
polynomial  tapering  function  was  employed  to  drive  the  induced-dipole/induced-dipole 
interactions  to  zero  at  this  cutoff,  with  scaling  starting  at  10.2  A. 

A  combined  MD-Monte  Carlo  (MD-MC)  approach30  31  was  utilized  to  obtain  crystal  relaxation 
at  1  atm  pressure  and  the  appropriate  temperatures  for  comparison  to  experimentally  available 
crystal  structures.  In  all  cases,  the  experimentally  derived  crystal  coordinates  were  used  as  a 
starting  configuration  for  the  MD-MC  runs.  The  MD-MC  technique  consists  of  two  iterating 
steps:  First,  a  500fs  trajectory  segment  of  canonical  ensemble  MD  simulation  (NVT-MD)  is 
performed,  followed  by  a  series  of  10  attempted  changes  to  the  simulation  cell  shape  and  volume 
via  isobaric-isothermal  rigid  molecule  Monte  Carlo  (NPT-MC)  moves  that  consisted  of  changing 
cell  parameters.  The  combined  NVT-MD  and  NPT-MC  steps  comprise  a  single  iteration,  with 
the  input  state  of  each  iteration  coming  from  the  results  of  the  previous  iteration.  A  total  of  5000 
iterations  were  performed  for  each  MD-MC  run,  generating  50000  attempted  MC  volume 
changes  and  2.5ns  of  composite  trajectory  time.  All  crystal  systems  were  simulated  for  at  least  2 
such  production  runs,  yielding  a  minimum  total  simulation  time  of  5  composite  nanoseconds.  In 
all  cases,  the  evolution  of  the  lattice  vectors  and  angles  were  monitored  over  time  to  ensure  that 
they  have  stabilized  at  representative  values  for  the  system  being  simulated,  with  simulation 
calculations  taken  only  from  the  latter  2ns  or  more  of  the  simulation  after  the  lattice  quantities  are 
stabilized. 
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IV.  Simulations  Results  and  Discussion 


Liquid  Density  and  Heat  of  V aporization  of  Neutral  Precursors.  Comparisons  between  measured 
values  and  MD  simulation  results  for  liquid  density  and  heat  of  vaporization  of  the  neutral 
hydrazines,  as  well  as  the  precursor  1-Azidobutane  used  to  parameterize  our  azide  anion,  are 
presented  in  Table  1.  All  of  the  molecules  have  been  investigated  at  298K  and  1  atm  of  pressure, 
with  the  resulting  predictions  for  both  density  and  Hvap  being  quite  accurate.  Particular  attention 
has  been  paid  to  accurate  reproduction  of  the  decreasing  trend  of  Hvap  with  the  addition  of 
subsequent  methyl  groups  in  the  hydrazine  system.  For  the  azide  group,  in  addition  to  fitting  well 
both  the  density  and  Hvap  of  1-Azidobutane,  the  azide  group  in  anionic  form  has  also  recently 
been  found  to  accurately  reproduce  densities  and  dynamic  properties  in  concert  with 
inridazoliunr-based  cations.32 

We  have  also  calculated  the  density  response  for  the  neutral  hydrazine  molecules  as  a  function 
of  temperature,  which  are  shown  in  Tables  2.  For  hydrazine,  density  parity  at  room  temperature 
is  quite  good,  but  as  the  system  heats  the  simulated  systems  expand  more  quickly  than 
experimental  observations.  This  effect  remains  modest  in  the  liquid  regime,  but  becomes 
significant  above  the  boiling  point  (387  K)  of  Hydrazine.  For  the  singly  methylated  molecule, 
we  obtain  reasonable  density  agreement  (-2.5%)  for  the  ranges  tested,  which  all  reside  below  the 
boiling  point  of  the  liquid  (361  K).  Even  better  agreement  between  simulation  and  experiment  is 
seen  for  UDMH,  with  simulations  up  to  slightly  above  the  boiling  point  of  336  K  falling  within  a 
density  variation  of  -1.5%.  The  trend  of  increasing  accuracy  with  increasing  nrethylation  of  the 
hydrazine,  coupled  with  the  extensive  fitting  which  has  been  performed  for  the  methyl  group 
employed,  indicates  a  small  deviation  in  the  force  field  fitting  for  the  NH2  group  of  hydrazine, 
which  is  very  quickly  and  easily  mitigated  by  the  presence  of  the  relatively  small  and  apolar 
methyl  group.  This  in  turn  implies  that  in  situations  where  the  local  environment  of  the  NH2 
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group  is  less  than  or  similarly  polarizable  as  that  seen  for  pure  hydrazine,  we  would  expect  the 
repulsion-dispersion  level  description  of  the  compounds  to  provide  an  adequate  definition  for 
good  reproduction  of  NH2  group  interactions. 

Crystal  Lattice  Parameters.  We  have  simulated  seven  test  crystals:  [Hyd][N03],33 
[MeHyd][N03],34  [DiMeHyd][N03],34  [MeHyd][N3],16  [DiMeHyd][N3],17  [Hyd][N(CN)2],  and 
[Hyd][CN7].15  For  the  entire  test  set  simulated,  we  find  generally  good  agreement  between  the 
simulated  crystals  and  those  obtained  from  experimental  investigations,  as  summarized  in  Table 
3.  With  the  exceptions  of  [Hyd][N03]  and  [MeHyd][N3],  cell  volumes  vary  by  less  than  3%,  with 
the  former  deviating  by  -6%  and  -7.1%,  respectively.  Individual  lattice  vectors  obtained  from 
MD  simulations  show  deviations  from  experiments  generally  falling  in  the  range  of  6%  absolute 
deviation  or  less,  while  lattice  angles  universally  deviate  by  a  roughly  2%  absolute  magnitude  or 
less.  Interestingly,  removing  an  extended  charge  center  (lone  pair)  on  the  N2H5+  cation  improves 
the  reproduction  of  the  [Hyd][N03]  crystal  structure.  However,  removing  the  extended  charge 
from  the  description  of  the  (CH3)N2H4+  cation  results  in  an  unstable  [MeHyd][N03]  crystal,  as 
shown  in  Figure  4.  One  of  the  important  features  of  the  developed  force  field  is  its  transferability 
and  subsequent  applicability  over  a  large  range  of  systems.  Therefore,  we  chose  to  systematically 
add  the  extended  charge  to  all  the  hydrazium-based  cations  because  it  significantly  improves  the 
electrostatic  potential  description  despite  obtaining  a  slightly  worse  description  of  the 
[Hyd][N03]  crystal  compared  to  the  one  for  the  model  without  the  extended  charge  on  N2H5+.  In 
general  we  have  found  that  the  presence  of  explicit  lone  pairs  (extended  charges)  is  necessary  in  a 
number  of  the  crystals  to  properly  dictate  the  location  of  the  NH2  group  hydrogens,  which  in  turn 
seems  to  be  necessary  for  the  formation  of  hydrogen  bonds  which  are  a  primary  component  of 
maintaining  crystal  cohesion  within  many  of  these  systems. 
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For  the  [Hyd][N(CN)2],  the  structural  symmetry  of  the  crystal  presents  a  challenge  for  both  the 
simulation  of  the  system  as  well  as  the  refinement  of  the  X-Ray  data.  The  final  crystal  structure 
is  in  the  Cmcm  space  group,  which  predicts  a  symmetric  orientation  for  the  hydrazinium  cations 
with  a  50%  fractional  occupancy  of  the  excess  cationic  proton  on  each  nitrogen  of  the 
hydrazinium  ion.  This  is  an  indicator  that  the  relative  orientation  of  the  NH3+  end  of  the 
hydrazinium  is  effectively  randomized  in  the  crystal,  with  neither  direction  being  preferred  over 
the  other  in  a  perfect  crystal.  Since  fractional  occupancy  is  a  statistical  effect  due  to  symmetry 
considerations,  we  are  forced  when  simulating  to  settle  upon  some  finite  implementation  of  the 
random  hydrazinium  orientations. 

Earlier  calculations  of  the  [Hyd][N(CN)2]  crystal  structure  indicated  a  Pnma  structure,  whose 
lattice  differs  from  that  of  the  Cmcm  structure  primarily  in  the  prediction  a  slightly  larger  tilt  of 
the  N-N  hydrazinium  axis  within  the  crystal  cell,  yielding  very  similar  lattice  vector  magnitudes 
and  unit  cell  volume.  It  is  from  this  data  which  we  took  the  stalling  point  for  our  simulations. 

The  presence  of  the  skew  within  the  Pnma  starting  point,  coupled  with  the  fact  that  the  proton 
location  are  primarily  predicted  from  electron  density  maps  and  symmetry  considerations,  leads 
to  unambiguous  placement  of  the  protons  within  the  Pnma  symmetry,  compared  to  the  fractional 
occupancy  of  the  Cmcm  group.  In  practical  terms,  simulation  from  the  Pnma  space  group  starting 
point  reflects  a  systematic  choice  of  excess  proton  placement  coupled  with  a  slight  additional 
skew  relative  to  the  correct  Cmcm  starting  point. 

For  experimental  comparison,  we  present  lattice  parameters  from  both  the  initially  calculated 
Pnma  crystal  structure,  as  well  as  the  final  calculation  yielding  a  Cmcm  crystal,  as  well  as  the 
associated  delta’s  for  the  simulation  with  respect  to  each  crystal.  In  the  case  of  calculating  the 
delta’s  for  the  Cmcm  crystal,  the  as-provided  lattice  vectors  represented  a  permuted  set  of  lattice 
vectors  with  respect  to  those  simulated  (and  calculated)  for  the  Pnma  crystal.  For  purposes  of 
ease  of  comparison,  we  have  removed  this  permutation  by  transforming  the  Cmcm  axes  as 
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follows.  ^fl^Cmcm^^Pnma)  </}>Cmci]i  <^  '>Pr!]'a-  tUld  <C > Cmcra^^^^Pnma-  Tills  pCllllUttition  WUS 
verified  by  visual  comparison  of  the  structure  file  of  the  Pnma  crystal  with  that  of  the  Cmcm 
crystal,  and  is  wholly  unambiguous.  We  compare  the  simulation  with  the  properly  permuted 
lattice  vectors  of  the  Cmcm  crystal  structure  in  Table  3.  Despite  the  slight  deviation  in  starting 
conditions,  it  can  be  seen  that  excellent  reproduction  of  the  crystal  structure  is  maintained, 
indicating  that  the  as-simulated  Pnma  crystal  is  indeed  only  a  slightly  perturbed  Cmcm  structure 
in  the  simulated  system  as  well  as  in  the  experimental  system. 

Structural  Comparisons  We  have  investigated  the  effect  of  increasing  methylation  on  the  ion 
correlation  structure  for  the  series  of  crystals  utilizing  the  NO3  anion  by  calculating  the  relative 
density  amplification  of  the  oxygen  of  the  N03  ion  near  the  hydrazinium  cations  as: 

p'  =  p{x,y,z)/pbulk  (12) 

where  p[x,y,z)  represents  the  local  density  of  the  target  atom  in  a  cubic  sub-volume  of 

approximately  0.1  A  per  side.  For  the  hydrazinium  cations,  we  define  the  local  coordinate  system 
relative  to  a  basis  formed  by  the  position  of  the  two  nitrogens  and  the  “outward”  lone  pair  force 
center  constant  (where  the  “outward”  lone  pair  is  that  which  simultaneously  has  a  larger  than  90° 
angle  with  each  of  the  NH2  N — H  bonds  as  well  as  the  N — N  bond).  For  notational  convenience, 
it  will  also  be  of  use  to  define  the  “top”  of  the  molecule  as  the  half  of  space  which  contains  the 
aforementioned  “outward”  pair  and  is  bisected  by  a  plane  embedding  the  N — N  bond,  with 
normal  perpendicular  to  the  N — N  bond  and  residing  in  the  plane  that  embeds  both  the  N — N  and 
N — “outward”  Lp  bond. 

We  begin  by  investigating  the  differences  in  the  hydrogen  bonding  structure  seen  in  the 
crystalline  systems  comprised  of  the  N2H5+  cation  with  all  of  the  available  anions  for  the 
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hydrazinium:  NO3',  CN7",  and  N(CN)2".  Figure  5  presents  graphs  of  the  3D  density  isosurfaces 
defined  in  Eqn.  12  with  p'  =  10,  with  the  differing  surface  colors  representing  the  negatively 
charged  atomic  species  for  their  associated  anions  (red  =  O  in  N03\  blue  =  N  in  N(CN)T,  yellow 
=  N  in  CN7 ),  with  differing  brightnesses  qualitatively  denoting  the  relative  charge  magnitude 
(brighter  is  more  negative).  Thus,  there  is  a  single  red  for  all  the  O  in  NO)  since  they  are  charge 
equivalent  (q  ~  -0.66e),  the  ligher  blue  represents  the  terminal  nitrogens  in  N(CN)T  (q  ~  -0.68e) 
while  the  dark  blue  represents  the  central  nitrogen  (q  ~  -0.62e),  and  there  are  three  shades  of  gold 
representing  the  terminal  azide  nitrogen  (bright,  q  ~  -0.36e),  the  bonded  azide  nitrogen  (medium, 
q  ~  -0.28e),  and  the  2,3  ring  nitrogens  (dark,  q  ~  -0.25e).  We  first  looked  for  areas  demonstrating 
similar  hydrogen  bonding  across  all  three  species,  rationalizing  that  these  would  represent  the 
primary  binding  sites  which  drive  the  crystalline  interactions.  In  order  to  more  clearly  see  just 
these  specific  locations,  Figure  5a  presents  the  density  isosurfaces  for  just  the  shared  isosurface 
locations,  with  all  additional  density  isosurfaces  removed.  The  overwhelming  importance  of 
strong  hydrogen  bonding  is  demonstrated  in  Figure  5a,  as  all  but  one  of  the  sites  highlighted  by 
the  requirement  of  concommittant  bonding  across  all  three  anion  species  directly  relate  to 
interactions  between  the  anions  and  the  hydrogens  of  the  N2H5+  molecule.  It  should  be  noted  that 
the  depiction  of  the  centralized  N2H5+  ion  is  provided  to  help  anchor  the  reference  of  the 
surrounding  density  isosurfaces.  In  reality,  the  NH3  group  including  the  positively  charged 
nitrogen  is  at  least  partially  mobile  (leading  to  the  production  of  charge  binding  surfaces  rather 
than  points )  for  the  N03  and  CN7  anions,  while  simulation  data  suggests  that  the 
[N2H5][N(CN)2]  crystal  the  NH3  group  is  able  to  fully  sample  it’s  conformational  rotation  in  the 
simulated  crystal.  Whether  this  is  a  mis-prediction  of  the  force-field  or  represents  an  increased 
inherent  mobility  of  this  group  due  to  intrinsic  differences  with  the  interaction  of  the  hydrogens 
and  N(CN)2‘  when  compared  to  the  other  anions  is  not  clearly  understood. 
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Figures  5b-d  present  a  number  of  different  viewpoints  on  the  anion  binding  with  the  same 
value  of  p\  with  the  principal  binding  sites  of  Figure  5a  left  solidly  colored,  while  ancillary 
binding  sites  arc  left  transparent.  In  this  case,  we  believe  all  of  the  transparent  sites  (some  of 
which  show  correspondence  between  differing  anions,  but  none  of  which  present  correspondence 
in  a  singularly  predictable  fashion  between  any  two  anions)  represent  the  effective  binding  due  to 
steric  and  packing  constraints.  These  sites  represent  the  inevitable  minimization  of  energy 
subject  to  the  constraints  of  intramolecular  connectivity  and  steric  intermolecular  forces  in 
conjunction  with  the  electrostatic  and  repulsion/dispersion  interactions  between  nearby  species. 

While  it  is  difficult  to  tell  just  from  this  type  of  analysis  what  the  overall  relative  contributions 
of  any  specific  interactions  are,  we  believe  the  principle  interactions  highlighted  in  Xa  and  arising 
from  strong  hydrogen  bonding  are  the  most  fundamental  in  determining  the  crystal  structure.  As 
further  proof  of  this,  we  consider  the  cases  of  both  the  (CH3)N2H4+  and  (CH3)2N2H3+  cations, 
where  we  have  two  crystals  for  each  with  the  nitrate  and  azide  counter-ions.  In  neither  case  are 
we  able  to  provide  graphs  such  as  those  shown  in  Fig.  5,  because  for  each  of  these  crystals  one  of 
the  two  crystalline  forms  stabilizes  in  a  manner  which  includes  a  conformational  state  which 
perturbs  the  torsional  state  of  the  methylated  nitrogen  from  it’s  lowest  energy  state 
([(CH3)N2H4][N3]  and  [(CH3)2N2H3][N03]  for  the  methylated  and  dimethylated  cation, 
respectively).  The  freezing  in  of  a  higher  conformational  energy  state  to  allow  more  facile 
hydrogen  bonding  is  indicative  of  the  overall  importance  of  the  hydrogen  bonding  interaction  in 
stabilizing  the  crystalline  system.  Given  the  importance  of  the  hydrogen  bonding  in  these 
systems,  we  have  also  investigated  the  relative  effect  that  methylation  has  upon  the  hydrogen 
bonding  by  looking  at  the  structural  correlations  in  the  nitrate  containing  crystals. 

To  understand  the  effect  of  progressive  methylation,  we  investigate  correlation  isosurfaces  for 
the  nitrate  bearing  crystals  by  examining  relative  density  amplification  of  the  nitrate  oxygens  at  a 
level  20  times  that  of  the  bulk  concentration.  In  the  [N2H5][N03]  system  shown  in  Fig.  6a,  there 
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is  ample  evidence  of  strong  association  between  the  oxygen  and  hydrogens,  with  the  positions 
near  the  triply  hydrogenated  end  of  the  N2H5+  cation  coordinating  an  oxygen  at  the  face  of  the 
pyramidal  base  formed  by  the  three  hydrogens,  as  well  as  showing  a  tendency  to  coordinate 
individually  to  these  hydrogens  in  the  roughly  triangular  series  of  “double”  density  regions 
located  commensurate  with  the  position  of  these  highly  charged  hydrogens.  Of  particular  interest 
is  the  fact  that  the  bottom  portion  of  these  “doubled”  regions  actually  shift  toward  the  NH2  end  of 
the  molecule  as  they  descend  to  take  advantage  of  the  extra  coordination  provided  by  the 
hydrogens  on  the  NH2,  despite  the  fact  that  the  NH2  hydrogens  also  coordinate  with  another  area 
directly  below.  Finally,  on  the  upper  portion  of  the  NH2  end  of  the  ion  a  somewhat  complex 
coordination  pattern  emerges,  presumably  due  to  the  necessity  to  gain  some  coordination  energy 
above  the  nitrogen  while  avoiding  the  electronegative  “well”  discussed  previously. 

When  the  N2H5+  cation  is  methylated  to  form  the  (CH3)N2H4+  cation,  the  resultant 
coordination  becomes  more  centralized  due  to  the  steric  interference  brought  about  by  the 
presence  of  the  methyl  group  (Fig.  6b).  In  this  and  the  following  (CH3)2N2H3+  (Fig.  6c)  graphic, 
it’s  important  to  remember  that  the  N+  group  has  a  symmetric  isomer  mirrored  across  the  plane 
which  includes  the  N-N  bond  and  bisects  the  HNH  angle  of  the  NH2  group,  where  in  both  cases 
the  bottom  substituted  group  is  constant  while  the  top  substituted  group  has  a  symmetric  parity 
across  this  mirror  plane.  Keeping  this  in  mind,  then,  it  is  easy  to  see  that  coordination  in  the 
[(CH3)N2H4][N03]  system  occurs  with  the  oxygen  interacting  with  two  hydrogens  from  the 
methyl  group  and  one  from  the  N+  substituted  hydrogen  on  the  upper  half,  or  the  two  hydrogens 
from  the  N+  substituted  hydrogen  and  one  from  the  gauche  hydrogen  on  the  NH2  group.  The  only 
additional  area  of  high  density  is  near  the  NH2  group,  situated  below  the  explicit  lone  pair  and 
above  the  NH2  hydrogens,  with  the  same  symmetry  as  the  groups  near  the  triply  substituted 
nitrogen. 
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We  have  also  investigated  the  difference  between  the  coordination  of  the  liquid  and  crystalline 
systems  where  possible.  For  [N2H5][CN7]  decomposition  occurs  within  the  solid  phase  before 
melting  occurs,  while  for  the  [N2H5][N(CN)2]  system  facile  isomerization  to  guanazole  occurs  at 
temperatures  well  below  the  melting  temperature,  leading  to  the  expectation  that  pure  [N2H5][ 
N(CN)2]  liquid  is  nearly  or  totally  non-existent.35  For  [N2H5][N03],  while  a  liquid  phase  is 
known  to  exist  experimentally,  we  were  unable  to  simulate  such  a  system  due  to  overly  strong 
interactions  between  the  acidic  hydrogen  atoms  attached  to  the  N+  atom  and  either  the  lone  pairs 
on  the  other  end  of  a  neighboring  hydrazinium  atom,  or  between  the  oxygen  of  the  [NO3]  ion  and 
the  acidic  hydrogen  attached  to  the  N+  atom.  For  the  remaining  systems  involving  the 
(CH3)N2H4+  and  (CH3)2N2H3+  cations  in  concert  with  the  N03  and  N3"  anions,  we  compare  the 
structure  from  the  crystal  at  experimentally  observed  temperatures  with  that  of  a  liquid  which  is 
uniformly  held  at  353K  across  all  liquid  samples.  The  liquid  simulation  was  run  for  a  minimum 
of  6ns,  and  the  local  density  enhancement  of  ion  correlations  calculated  as  above. 

Figure  7  presents  a  comparison  of  the  observed  ion  correlation  structures  for  the  liquid  and 
crystalline  systems  for  the  four  systems  for  which  liquid  data  was  obtainable.  For  both 
[(CH3)N2H4][N03]  and  [(CH3)2N2H3][N3]  (Figs.  7a  and  7d)  the  differences  between  crystal  and 
liquid  packing  are  rather  small,  with  the  relaxed  steric  constraints  of  the  liquid  phase  allowing  the 
anions  to  interact  primarily  through  the  acidic  hydrogens  attached  directly  to  the  respective 
nitrogens  in  the  liquid,  rather  than  forcing  interaction  with  the  only  slightly  charged  hydrogens  of 
the  methyl  groups,  as  occurs  in  the  crystal  phase.  By  comparison,  both  [(CH3)2N2H3][N03]  and 
[(CH3)N2H4][N3]  (Figs.  7b  and  7c)  show  significantly  larger  differences  between  the  crystal  and 
liquid  phases.  For  both  of  these  cations,  in  the  crystalline  phase  there  is  sufficient  binding  energy 
due  to  the  proximity  of  the  anions,  and  the  hydrogen  bonding  which  accompanies  this  proximity, 
to  lead  to  the  stabilization  of  cationic  conformations  in  the  crystal  phase  which  are  not  the  lowest 
energy  conformer.  In  both  cases,  the  presence  of  external  hydrogen  bonds  leads  to  a 
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conformation  in  the  crystal  phase  with  a  methyl  group  in  the  trans  position  relative  to  the  upward 
lone  pair  of  the  NH2  group  (represented  by  the  transparently  depicted  groups  in  Figs.  7b  and  7c), 
which  serves  the  effect  of  exposing  the  acidic  hydrogens  on  the  substituted  N+  atom  to  attack 
from  anionic  groups  from  all  angles,  allowing  the  crystal  to  stabilize.  In  the  liquid  phase,  where 
accommodation  of  strong  steric  interactions  loses  importance,  the  cations  regain  their  standard 
low  energy  conformation  with  hydrogens,  rather  than  methyl  groups,  in  the  trans  position.  This  is 
especially  notable  in  the  [(CH3)N2H4][N3],  where,  keeping  the  N+-N-Lp  basis  as  constant,  the 
difference  in  conformation  leads  to  a  near  total  inversion  in  the  spatial  location  of  primary  anion 
association  near  the  N+  atom.  For  both  of  these  systems,  the  net  effect  is  a  far  larger  localization 
of  anionic  counter-charge  in  the  liquid  phase  compared  to  the  crystalline  phase. 

Mechanical  Properties 

We  utilize  the  formalism  of  Parinello  and  Rahman,36  which  relates  the  elastic  stiffness 

C  (b-£  ) 

tensor  (  ljkl)  to  the  fluctuations  in  the  elastic  strain  tensor  \  lJ  u  >  of  the  material,  to  calculate  the 

elastic  behavior  of  our  simulated  crystals  via  the  observed  fluctuations  in  the  lattice  vectors  from 

our  MD/MC  simulations.  In  the  equations  that  follow,  we  use  a  bolded,  double-overbarred 

quantity  ( A )  to  represent  a  four-dimensional  tensor  and  a  bolded,  single-overbarred  quantity  ( A 
)  to  represent  a  two-dimensional  tensor  (matrix).  A  bolded  quantity  ( A )  is  a  vector,  and  a  non- 

bolded  quantity  ( A)  a  scalar.  The  inner  product  is  represented  by  A  :  B  while  an  outer  product 
is  represented  by  A®  B . 

In  a  generic  elastic  solid,  the  fundamental  stress-strain  equations  are:37 


o  =  C  :e 

(13) 

£  =  S  :o 

(14) 
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where  cr  is  the  stress  tensor,  £  the  strain  tensor,  S  the  compliance  tensor,  and  C  the  stiffness 
tensor. 

Due  to  the  symmetry  of  the  stress  and  strain  tensors,  it  is  possible  to  rewrite  Eqns.  13  and 
14  in  a  contracted  form,  utilizing  Voigt  notation  ( xu  =  x{,  x12  =  x2,  x33  =  x3,  x23  =  x32  =  x4, 
xl3  =  x3l  =  x4,  x12  =  X21  =  x5 ).  This  contracts  cr  and  £  each  into  a  6-dimensional  vector,  while 

S  and  C  become  the  familiar  6x6  matrices  normally  encountered  when  investigating  elastic 
constants.  Due  to  ambiguity  in  the  possible  definitions  of  the  contracted  stress  and  strain,  the 
contracted  stiffness  and  compliance  matrices  are  not  necessarily  matrix  inverses  of  one  another. 
However,  if  one  uses  the  construction  of  Tsai,38  and  defines: 

cr  —  (<Jll,o22,033,023,cri3,crl2)  (15) 

the  resulting  contraction  of  the  compliance  and  stiffness  tensors  yield: 

SvtJvkl  =  Sm  for  {l  <  Vi;  <  3  and  1  <  Vkl  <  3} 

SVlJvtl  =  2 Sijki  for  {l  <  Vu  <  3  and  4  <  Vkl  <  6}  or  [4  <  Vtj  <  6  and  1  <  Vkl  <  3}  (17) 

SvtJvkl  =  4  Sijki  for  {4  <Vij<6  and  4  <  Vkl  <  6} 

CvijVu  =  Cm  for  all  VipVkl  (18) 

Under  these  definitions,  the  contracted  matrices  C  and  S  are  simple  matrix  inverses  of  each 
other. 

Calculation  of  elastic  constants  may  then  proceed  by  noting  that:36 

S  =  [e®~e)^  (19) 

where  the  angle  brackets  denote  system  averages,  S  is  the  compliance,  V  is  the  volume  of  the 
simulated  cell  at  temperature  T,  ^ b  is  Boltzmann's  constant,  and  £  is  the  strain,  given  by: 
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e  =  -§(h'oW-l) 


(20) 


I_  - - 

with  no  is  the  reference  state  of  the  system,  and  G  =  h'h  is  the  metric  tensor  of  the  current 

(instantaneous)  system  state.  Here,  a  prime  (  '  )  represents  a  matrix  transpose.  We  take  no  to  be 
equivalent  to  the  average  over  the  entire  set  of  sampled  cells. 

We  have  performed  these  calculations  on  all  7  of  our  simulated  crystals,  with  the  results 
summarized  in  Table  4  for  the  principal  compression  ( C, , ,  C22 ,  C33 )  and  shear  ( C44 ,  C55 ,  C66 ) 

coefficients.  Since  simulations  occurred  at  different  temperatures,  a  direct  comparison  amongst 
the  entirety  of  the  data  is  somewhat  problematic.  The  two  significant  outliers  for  temperature  are 
the  [N2H5][N03]  system  at  120K  and  the  [N2H5][CN7]  system  at  100K.  The  remainder  of  the 
systems  fall  at  either  173K  or  200K,  both  of  which  are  close  enough  relative  to  the  other  to  expect 
qualitative  trends  to  be  accurate.  In  addition,  it  is  reasonable  to  expect  that  the  experimental 
conditions  for  the  unambiguous  determination  of  crystal  structure  represents  a  rough  regime  of 
“normal”  solid  handling  conditions  which  might  be  representative  of  the  behavior  of  the 
crystalline  systems  in  actual  utilization.  Under  the  second  assumption,  we  can  therefore 
qualitatively  compare  the  stiffness  response  of  the  crystals  in  order  to  look  for  broad  general 
trends  in  mechanical  behavior.  For  the  compression  stiffness,  we  find  [N2H5][CN7]  ~ 
[N2H5][N(CN)2]  >  [N2H5][N03]  >  [(CH3)N2H4][N03]  ~  [(CH3)N2H4][N3]  >  [(CH3)2N2H3][N03]  ~ 
[(CH3)2N2H3]  [N3] . 

For  shear  stiffness  the  relations  are  more  complex.  In  all  cases  the  shear  stiffness  is  smaller 
than  the  compression  stiffness,  with  relative  compressive  to  shear  stiffness  ratios  roughly  ranging 
from  2  to  10,  depending  on  the  system  and  direction.  The  size  of  the  anion  and  its  ability  to  form 
a  hydrogen  bonding  network  seems  to  correlate  well  with  the  shear  stiffness  of  ionic  crystals.  For 
example,  [(CH3)2N2H3][N3]  with  the  smallest  anion  and  the  least  number  of  hydrogen  bonds  per 
ion  pair  (or  per  volume)  has  the  lowest  stiffness,  replacing  a  methyl  group  with  the  hydrogen 
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((CH3)2N2H3+  — >  (CH3)N2H4+)  increases  the  number  of  hydrogen  bonds  per  ion  pair  and  increases 
the  stiffness.  This  trend  is  clearly  observed  for  the  [N2H5][N03]  >  [(CH3)N2H4][N03]  > 
[(CH3)2N2H5]  [N03]  ionic  crystals.  Increasing  the  size  of  anion  from  N3  to  N03 ,  N(CN)2"  and 
CN7  for  the  (CH3)N2H4+-based  ionic  crystals  generally  results  in  a  larger  shear  stiffness  for  the 
larger  anions,  with  the  most  anisotropic  anion  CN7‘  yielding  the  most  anisotropic  shear  stiffness 
matrix.  When  comparing  the  crystals  involving  [N03]  and  [N3],  it  is  interesting  to  note  that  while 
both  anions  produce  similar  compression  strengths,  there  is  a  substantial  difference  between  their 
shear  strengths,  perhaps  due  to  the  different  shape  of  the  planar  (N03)  group  and  the  linear  (N3)\ 

V.  Conclusions 

Utilizing  a  many-body  polarizable  formulation  with  explicit  lone  pairs,  we  have  parameterized 
a  transferrable,  consistent  force  field  to  represent  the  hydrazinium  cation  and  its  mono-methyl  and 
unsymmetrically  dimethylated  derivatives  in  conjunction  with  N3",  N03\  N(CN)2"  and  CN7~  anion. 
We  have  simulated  seven  hydrazinium-based  crystals  across  four  differing  anionic  species,  with 
most  predictions  showing  excellent  volumetric  parity  with  experiment,  on  a  range  of  -2.5%  or 
less,  with  individual  lattice  vector  deviation  of  -1-6%.  For  the  [N2H5][N03]  crystal  volumetric 
deviation  is  -6%,  and  for  the  [(CH3)N2H4][N3]  crystal  -7%,  with  individual  lattice  vectors  in  the 
1-7%  range.  Of  particular-  interest  in  the  case  of  the  cations  is  the  necessity  of  including  lone 
pairs  to  accurately  reproduce  the  local  electrostatic  field  of  these  ions  near  the  NH2  group.  While 
lone  pair  inclusion  is  relatively  standard  for  both  anions  and  certain  neutral  molecules,  the 
necessity  of  their  presence  on  a  cationic  species  is  noteworthy.  We  have  also  performed  liquid 
state  simulations  of  the  neutral  hydrazines,  as  well  as  the  [N03]  and  [N3]  based  ionic  salts. 
Comparison  of  the  hydrazines  shows  excellent  reproduction  of  both  density  and  heat  of 
vaporization  at  room  temperature,  with  density  reproduction  drifting  for  the  pure  hydrazine  at 
elevated  temperatures,  while  predictions  for  [MH]  and  [UDMH]  remain  stable. 
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In  ion  pair  QC  calculations  at  the  MP2//cc-pvTz  theory  level,  [N2H5HNO3],  [N2H5][CN7], 
[(CH3)N2H4][N03],  [(CH3)N2H4][N3],  and  [(CH3)2N2H4]  [N3]  all  show  an  energetic  ground  state 
for  each  pair  with  the  (cation)  anion  (de)protonated  in  the  isolated  gas  phase.  For  the 
[N2H5][N03]  pair  calculation,  we  have  also  done  QC  in  a  polarizable  continuum,  noting  that 
placing  the  de-ionized  pair  into  a  polarizable  continuum  leads  to  re-ionization  of  the  constituent 
neutral  molecules.  We  believe  this  is  indicative  of  the  fact  that  while  isolated  ion  pairs  are 
unstable,  bulk  phase  ionic  liquids  are  effectively  stabilized  in  the  ionic  state  due  to  the  dielectric 
medium  of  the  surrounding  ionic  liquid. 

It  is  noteworthy  that  of  the  seven  crystals  we  simulate,  we  see  the  largest  deviation  in  those 
that  have  the  smallest  anions.  Furthermore,  the  quality  of  prediction  seems  to  systematically 
increase  with  increasing  methylation  of  the  hydrazinium  cations.  Coupled  with  the  above 
discussions  of  results  in  the  QC  calculations,  and  noting  that  we  parameterize  the  charge  and 
polarizability  each  from  the  isolated  gas  phase,  we  believe  it  may  be  likely  that  overly  strong 
charge  and/or  polarization  interactions  with  respect  to  the  acidic  hydrogens  on  the  N+  atoms  are 
responsible  for  the  undesirable  force  field  behavior.  As  the  molecular  structures  inherently  begin 
to  further  screen  these  interactions  (either  through  increased  methylation  of  the  cation  or 
increasing  anion  size,  both  of  which  impedes  close  approach  to  the  acidic  hydrogens  in  the  bulk), 
the  undesirable  behavior  goes  away.  This  would  be  consistent  with  the  observation  that 
[(CH3)N2H4][N03]  yields  better  results  than  [(CH3)N2H4][N3],  since  the  N3“  anion  is  smaller  than 
the  NO;/  anion  and  thus  more  able  to  attack  the  hydrogens  near  the  methyl  group  of  (CH3)N2H4+, 
as  well  as  the  better  quality  results  of  [(CH3)2N2H3][N3]  over  [(CH3)N2H4][N3]  and 
[(CH3)N2H4][N03]  over  [N2H5][N03]. 

Investigation  of  anion  bonding  similarities  with  a  common  N2H5+  across  three  differing 
crystals  demonstrates  the  importance  of  hydrogen  bonding  in  stabilizing  the  various  crystals.  We 
found  significant  differences  between  hydrazinium  cation  coordination  by  anions  in  the  crystal 
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and  liquid  phase  indicating  that  the  liquid  structure  should  be  inferred  from  the  crystal  structure 
with  caution.  Structural  features  such  as  the  hydrogen  bonding  network  and  ion  sizes  and  shapes 
were  correlated  with  the  shear  stiffness  matrix  of  ionic  crystals  calculations  from  MD 
simulations. 
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Table  1:  Comparison  between  experimental  and  simulated  densities  and  enthalpies  of 
vaporization  for  the  neutral  hydrazine  liquids. 


N2H4  (CH3)N2H3  (CH3)2N2H2  1-Azidobutane 

Temp.  (K)  298  298  298  298 

Dens,  (kg/nr3) 


Experiment 

Simulation 

1003.6 

1008.7 

870 

871 

791 

800 

oo  oo 
oo  oo 

i— *  OJ 

A(%) 

0.5 

0.1 

1.1 

0.2 

AHVap  (kcal/Mol) 

Experiment 

10.68 

9.68 

8.42 

9.17 

Simulation 

10.68 

9.77 

8.09 

9.27 

A(%) 

0.0 

0.9 

-3.9 

1.1 
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Table  2:  Comparison  between  predicted  and  observed  densities  of  the  hydrazine  liquid  force  field 


at  elevated  temperatures. 


Hydrazine 

T  (K)  299.75  366.47  449.81 


Experiment 

1003 

940 

858 

Sim. 

1007 

916 

111 

A(%) 

0.4 

-2.5 

-9.4 

Methylhydrazine 

T  (K) 

308.15 

333.15 

343.15 

Experiment 

858 

834 

824 

Sim. 

846 

815 

801 

A(%) 

-1.5 

-2.3 

-2.8 

N,N -Dimethylhydrazine 

T  (K) 

273.2 

310.26 

338.73 

Experiment  | 

812 

773 

742 

Sim. 

827 

786 

753 

A(%) 

1.8 

1.7 

1.5 
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Table  3:  Comparison  between  experimentally  observed  and  simulated  crystal  structures.  Lattice 
vectors  are  in  A,  lattice  angles  in  degrees,  and  cell  volume  in  A3. 
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<a>  <b>  <c>  <a>  <b>  <g>  <Vol.> 


[N2H5][N03] 


Experiment 

7.965 

5.657 

8.122 

90.00 

91.34 

90.00 

365.854 

Simulation 

7.44 

5.97 

7.75 

91.33 

92.55 

89.88 

343.370 

A(%) 

-6.6 

5.4 

-4.6 

1.5 

1.3 

-0.1 

-6.1 

[(CH3)N2H4][N03] 


Experiment 

3.779 

11.342 

11.107 

90.000 

99.086 

90.000 

470.147 

Simulation 

3.96 

10.71 

11.30 

90.08 

100.61 

90.04 

471.34 

A(%) 

4.9 

-5.6 

1.7 

0.1 

1.5 

0.0 

0.3 

[(CH3)2N2H3][N 

Experiment 

Simulation 

o3] 

14.039 

14.48 

5.649 

5.56 

7.603 

7.37 

90.000 

89.69 

90.000 

89.76 

90.000 

89.99 

602.958 

593.61 

A(%) 

3.1 

-1.5 

-3.1 

-0.3 

-0.3 

-0.0 

-1.6 

[N2Hs][CN7] 


Experiment 

10.811 

7.464 

7.668 

90.000 

101.437 

90.000 

606.469 

Simulation 

10.79 

7.47 

7.49 

89.78 

101.15 

90.08 

592.03 

A(%) 

-0.2 

0.0 

-2.3 

-0.3 

-0.3 

0.1 

-2.4 

[(CH3)N2H4][N3] 


Experiment 

9.963 

5.132 

9.228 

90.000 

90.000 

90.000 

471.773 

Simulation 

10.42 

4.84 

8.67 

89.88 

90.02 

90.03 

437.25 

A(%) 

4.6 

-5.6 

-6.1 

-0.1 

0.0 

0.0 

-7.3 

[(CH3)2N2H3][N3] 
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Experiment 

Simulation 

A(%) 

[n2h5][N(cn)2; 

Exp.  (Pnma) 
Exp.  (Cmcm) 
Simulation 

Apnma(%) 

Acmcm(%) 


8.445 

7.042 

10.075 

90.000 

108.321 

90.000 

568.741 

8.49 

7.18 

10.04 

88.50 

108.90 

90.69 

578.29 

0.5 

1.9 

-0.3 

-7.7 

0.5 

0.8 

7.7 

6.581 

9.000 

7.497 

90.000 

90.000 

90.000 

444.033 

6.568 

8.983 

7.484 

90.000 

90.000 

90.000 

441.559 

6.82 

8.89 

7.31 

90.06 

90.03 

89.88 

443.62 

3.7 

-7.2 

-2.5 

0.1 

0.0 

-0.1 

-0.1 

3.8 

-1.0 

-2.3 

0.1 

0.0 

-0.1 

0.5 
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Table  4.  Principal  elastic  constants  (in  GPa)  for  the  crystals  investigated  in  this  work. 


[N2H5] 

[(CH3)N2H4][(CH3)2N2H3] 

[N2H5] 

[(CH3)N2H3][(CH3)2N2H3] 

[N2H5] 

[no3] 

[no3] 

[N03] 

[cn7] 

[N3] 

[N3] 

[DCA] 

T 

120.0  K 

200.0  K 

200.0  K 

100.0  K 

200.0  K 

173.0  K 

173.0  K 

Cn 

42.70+1.21 

38.38+0.90 

19.20+0.60 

55.71+1.48 

44.83+0.71 

23.87+0.72 

54.00+1.31 

C22 

42.88+1.46 

38.95+0.92 

18.58+0.49 

97.73+2.82 

33.45+0.52 

15.03+0.54 

94.38+3.59 

C33 

109.12+3.09 

37.40+0.76 

11.60+0.37 

60.39+2.17 

22.98+0.49 

19.06+0.84 

43.95+1.77 

C44 

22.49+0.75 

21.88+0.25 

6.94+0.27 

6.19+0.17 

4.44+0.17 

2.54+0.07 

11.18+0.29 

C55 

21.13+0.64 

7.38+0.16 

6.82+0.20 

28.83+0.89 

4.94+0.09 

2.34+0.06 

12.56+0.52 

c66 

22.33+0.97 

9.62+0.21 

4.25+0.14 

9.90+0.17 

8.07+0.35 

3.64+0.13 

10.37+0.29 
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Figure  Captions 


Figure  1. 

Schematics  for  the  neutral  hydrazines  (a)  used  as  a  basis  for  the  hydrazinium  cations  (b)  as  well 
as  for  the  anions  (c)  employed  in  this  work. 

Figure  2. 

Representation  of  the  difference  in  the  shapes  of  the  force  field  derived  (solid)  and  quantum 
chemistry  calculated  (MP2/aug-cc-pvDz,  mesh)  electrostatic  surfaces  at  0.125eV  for  the 
hydrazinium  cation  without  (a)  and  with  (b)  explicit  lone  pairs.  The  complete  inability  of  the 
system  without  lone-pairs  to  even  qualitatively  capture  the  electrostatic  potential  well  near  the 
tetrahedral  position  of  the  NFF  group  is  evident,  with  the  force  field  derived  potential  surface  in 
part  (b)  clearly  descending  toward  the  lone  pairs  (represented  by  yellow  spheres). 


Figure  3: 

Comparison  of  potential  isosurfaces  at  a  value  of  -0.15eV  for  the  [CN7]  anion.  Quantum 
chemistry  calculations  are  represented  by  the  green  mesh,  while  those  for  the  forcefield  without 
(red  solid)  and  with  (light  blue,  solid)  lone  pairs  are  superimposed.  There  are  5  planar  lone  pair 
force  centers:  One  each  on  the  ring  nitrogens,  with  the  5th  attached  to  the  C-N-N  bend  attaching 
the  tetrazole  ring  to  the  azide  tail. 

Figure  4: 
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Comparison  of  individual  normalized  lattice  vector  evolution  for  the  [MeHyd][N03]  crystal 
system  for  simulations  which  include  (solid)  and  ignore  (dashed)  the  effect  of  explicit  lone  pair 
force  centers.  The  predicted  lattice  vector  normalized  by  the  experimental  value  is  plotted  for  the 
a(green),  b(blue),  and  c(red)  principle  lattice  directions. 

Figure  5: 

Ion  correlations  between  the  hydrazinium  cation  and  all  simulated  anions.  Density  isosurfaces 
represent  local  amplification  of  the  density  of  specific  atoms  of  the  anionsby  a  factor  of  10  over 
their  expected  bulk  values:  oxygen  in  N03  (red),  terminal  (light  blue)  and  central  (dark  blue) 
nitrogens  in  N(CN)2\  terminal  azide  (bright  yellow),  bonded  azide  (medium  gold),  and  N3/N4 
ring  (dark  gold)  nitrogens  in  CN7\ 

Part  (a)  presents  only  the  density  isosurfaces  for  locations  which  are  fundamentally  similar  for  all 
three  anions,  representing  the  principle  locations  of  hydrogen  bonding,  while  differing  views  are 
shown  by  (b  -  side),  (c  -  NH3  end)  and  (d  -  below)  with  the  surfaces  of  (a)  colored  opaquely 
while  secondary  locations  (all  at  the  same  relative  density)  which  differ  from  anion  to  anion  are 
represented  with  transparent  versions  of  the  appropriate  color. 

Figure  6: 

Ion  correlations  between  the  hydrazinium  cations  and  the  oxygen  in  the  nitrate  anions  for  (a) 
hydrazinium  nitrate,  (b)methylhydrazinium  nitrate,  and  (c)dimethylhydrazinium  nitrate.  Density 
isosurfaces  represent  local  amplification  of  oxygen  density  that  is  20  times  greater  than  the  bulk 
density  of  oxygen  in  the  crystal. 

Figure  7 : 
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Comparison  of  ion  correlations  for  crystalline  and  liquid  systems  for  (a)  [MeHyd][N03], 

(b) [DiMeHyd][N03],  (c)[MeHyd][N3],  and  (d)[DiMeHyd][N3].  Crystal  correlations  are 
represented  by  the  meshed  areas,  while  liquid  correlations  by  the  solid  yellow  areas.  For  (b)  and 

(c) ,  the  cations  undergo  conformational  changes  as  the  crystal  melts,  with  representative 
conformations  for  the  crystalline  system  indicated  by  transparency  in  the  functional  groups  which 
show  positional  deviation  due  to  the  differences  in  average  conformations. 

Figures 


Unsymmetrical  dimethylhydrazine 


Figure  la 
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Hydrazinium 

Methylhydrazinium 


Figure  lb 


Nitrate 


5-Azidotetrazolate 


Figure  lc 
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Figure  2b 
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Percent  Deviation  in  Lattice  Vector 


Figure  3 


MC  Attempt  (xIO3) 


Figure  4 
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Figure  5  a 


Figure  5b 


Figure  5  c 
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Figure  5d 


Figure  6a 
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Figure  6b 


Figure  6c 


Figure  7  a 
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Figure  7b 


Figure  7c 
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Figure  7d 
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Supplemental  Information 

Determination  of  [N^HdiNtCNEl  Crystal  Structure 

Experimental.  Hydrazinium  dicyanamide  was  prepared  according  to  literature  [1],  Single 
crystals  were  obtained  from  methanol  solution  layered  with  diethylether. 

X-ray  Analyses.  The  single-crystal  X-ray  diffraction  data  were  collected  on  a  Bruker  3- 
circle-platform  diffractometer  equipped  with  a  SMART  APEX  2  detector  with  the  y-axis  fixed  at 
54.74°  and  using  MoK„  radiation  from  a  fine-focus  tube.  The  goniometer  head,  equipped  with  a 
nylon  Cryoloop  and  magnetic  base,  was  used  to  mount  the  crystals  using  perfluoropolyether  oil. 
The  data  collection  as  well  as  structure  solution  and  refinement  were  canned  out  using  standard 
procedures  with  the  APEX2  V.2.1-4,  SMART  V.5.622,  SAINT  7.24A,  SADABS,  and  SHELXTL 
software  packages  and  programs. [2]  Crystal  data  and  refinement  details  of  crystals  of  hydrazinium 
dicyanamide  are  given  in  Table  1. 

[1]  Frankel,  M.B.;  Burns,  E.A.;  Butler,  J.C.;  Wilson,  E.R.  J.  Org.  Chem.  1963,  28,  2428. 

[2]  APEX2  V.2.1-4,  SMART  V.5.622,  SAINT  7.24A,  SADABS,  SHELXTL  ed.;  Bruker-AXS, 
INC.:  Madison,  WI  USA,  2007. 
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Table  SI  Crystal  and  structure  refinement  data  for  hydrazinium  dicyanamide. 


Hydrazinium  dicyanamide 

Formula 

c2h5n5 

Space  group 

Cmcm  orthorhombic 

a  (A) 

8.983(4) 

b  (A) 

7.484(3) 

c  (A) 

6.568(3) 

a,  (3,  y(°) 

90 

v/A3 

441.5(3) 

Pcalc./g  Cm'3 

1.491 

Z 

4 

Formula  weight 

99.11 

p/mm'1 

0.112 

Temperature  (K) 

173(2) 

X(MoKa) 

0.71073 

Crystal  size 

0.20x0.15x0.10 

Theta  range  0/° 

3.54  to  25.34 

Index  range 

-10<  h<10,  -8<k<8,  -7<1<7 

Reflection  collected 

2097 

Independent  [R(int)]/ 

234  [0.0170] 

Obs.  refl.  ([I  >  2.0  o(I)]) 

217 

F(OOO) 

208 

GooF 

1.186 

Ri,  w R  [I  >  2o(I)] 

0.0264,  0.0738 

Ri,  w R2  (all  data) 

0.0282,  0.0750 

L.diff.  peak/hole  eA3 

0.14  and  -0.17 

Absorption  correct. 

multiscan  SADABS 

T  T 

A  min?  A  max 

0.975,  0.990 

Data/restraints/par  am. 

234/0/30 

Refinement  method 

Full-matrix  least  squares  on  F2 

v r _ /it-'  |2  it-1  i2\2t  i  Vi 

CCDC-xxxxxx  contains  the  supplementary  crystallographic  data  for  this  paper.  These 
data  can  be  obtained  free  of  charge  from  the  Cambridge  Crystallographic  Data  Centre  via 
http :  //www.  cede .  cam .  ac .  uk/  data_reque  s  t/cif . 
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Figure  SI.  ORTEP  diagram  showing  conformation  and  atom  numbering  scheme  of 
hydrazinium  dicyanmide. 
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